function [Results] = pred_VARboot(y , X , c , lags , weight , Nboot , method , W , R2prc)
% This code runs multivariate regressions with bootstrap.
% We do the following
% 1. Estimate the regression with OLS and store the residuals
% 2. Estimate a VAR(1) for the predictors. Bias-adjust the coefficients and
% store the residuals (Bekaert and Hodrick 2001, Adam et al. 2017)
% 3. Simulate under the alternative and null
% 4. Re-run regressions to get coeff. bias, p-value and R2 distribution
%
%  INPUT:
%      y      - T by 1 vector
%      X      - T by K matrix, DO NOT LAG
%      c      - Whether to include a constant or not
%      lags   - Number of lags used in calculating s.d.
%      weight - s.d. matrix weighting. 1 for newey-west weighting 
%                                      0 for even weighting
%                                     -1 skips standard error computations.
%      Nboot  - Number of bootstrap simulations
%      method - Method of bootstrap data residuals
%             - 1: i.i.d. bootstrap
%             - 2: circular bootstrap
%             - 3: stationary bootstrap
%      trend  - Whether to include a trend in estimating VAR(1) for X[OPTIONAL]
%               Default is 0 not including a trend.
%      R2prc  - percentile of adj. R2 in the bootstrap null[OPTIONAL].
%               Default 95.
%  OUTPUT:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check inputs
if nargin < 7
    error('7 inputs required!')
end
% Get dimensions
[T_1 , K] = size(X); 
if size(y , 1) ~= T_1
    error('Dimensions not consistent!')
end
if size(y , 2) > 1
    error('Only univariate LHS allowed!')
end

T = T_1 - 1;%For predictive regressions
py = y(2 : end);%LHS to be predicted
LX = X(1 : end - 1 , :);%Lagged predictors

% Constant option
if c == 1
    Const = ones(T , 1);
elseif c == 0
    Const = [];
else
    error('Wrong constant choice!')
end

% Block length. Following Crump and Gospodinov (2019) to use max 
XL_optimal = opt_block_length(X); %See Andrew Patton's Website
if method == 1 % i.i.d
   XL = 1;
   type = 1;
elseif method == 2 % circular
   XL = ceil(max(XL_optimal(2 , :)));
   type = 2;
elseif method == 3 % stationary
   XL = ceil(max(XL_optimal(1 , :)));
   type = 3;
end

% Default options
switch nargin
   case 7
      R2prc = 95;
   case 8
      XL = W;
      R2prc = 95;
end

%% OLS estimates
RHS = [Const , LX];
[b_raw , se_raw , ~ , R2adj_raw] = olsgmm(py , RHS , lags , weight);

e_raw = py - RHS * b_raw; % Store residuals
dy = py - mean(py); % De-meaned y

t_raw = b_raw./se_raw; % Unadjusted t-statistics

sign_b = sign(b_raw);% For one-sided tests

%% Estimate bias-adjusted first order VAR for X
[~ , Xb_adj , Xe_adj] = VAR1_bsadj(X , Nboot , method , 0); 

%% Generate bootstrap samples
b_alt = nan(K + c , Nboot); % Store coefficients under the alternative
t_alt = nan(K + c , Nboot); % Store t under the alternative
t_null = nan(K + c , Nboot);% Store t under the null
R2adj_null = nan(1 , Nboot);% Store adj. R2 under the null

% Bootstrap residuals (index)
error_index = bs_choose((1 : T)' , Nboot , XL , type);%Map to 2 to T+1

% Bootstrap sample
X_temp = nan(T_1 , K);
y_alt_temp = nan(T_1 , 1);
y_null_temp = nan(T_1 , 1);
% Conditioning on the first observation (Kothari and Shanken 1997; Crump and Gospodinov 2019)
X_temp(1 , :) = X(1 , :); 
y_alt_temp(1) = y(1);
y_null_temp(1) = y(1);
% Each simulation
for i = 1 : Nboot   
    error_index_temp = error_index(: , i);%Map to 2 to T_1, bootstrapped
    Xe_temp = Xe_adj(error_index_temp , :);%Bootstrapped X shocks from 2 to T+1
    e_temp = e_raw(error_index_temp);%Bootstrapped y shocks from 2 to T+1
    % Bootstrap X          
    for j = 2 : T_1              
        X_temp(j , :) = [1 , X_temp(j - 1 , :)] * Xb_adj + Xe_temp(j - 1 , :);
    end
    RHS_temp = [Const , X_temp(1 : end - 1 , :)];%Lagged predictors
    % Simulate under the alternative
    y_alt_temp(2 : end) = RHS_temp * b_raw + e_temp; 
    % Simulate under the null
    y_null_temp(2 : end) = mean(y) + e_temp;           
    
    % Under the alternative the maintained null is b = b_raw
    [b_alt_temp , se_alt_temp] = olsgmm(y_alt_temp(2 : end) , RHS_temp , lags , weight);
    b_alt(: , i) = b_alt_temp;
    t_alt(: , i) = (b_alt_temp - b_raw)./se_alt_temp; 
        
    % Under the null the maintained null is b = 0
    [b_null_temp , se_null_temp , ~ , Radj_null_temp] = olsgmm(y_null_temp(2 : end) , RHS_temp , lags , weight);
    t_null(: , i) = (b_null_temp -[mean(y); zeros(K , 1)])./se_null_temp;
    R2adj_null(i) = Radj_null_temp;
end

% One-sided p-value
p1_alt   = mean(kron(t_raw .* sign_b , ones(1 , Nboot)) < t_alt .* repmat(sign_b , 1 , Nboot) , 2);
p1_null = mean(kron(t_raw .* sign_b , ones(1 , Nboot)) < t_null .* repmat(sign_b , 1 , Nboot) , 2);

% Two-sided p-value
p2_alt = mean(abs(kron(t_raw , ones(1 , Nboot))) < abs(t_alt) , 2);
p2_null = mean(abs(kron(t_raw , ones(1 , Nboot))) < abs(t_null) , 2);

% Adjusted coefficients
mb = mean(b_alt , 2); 
bias = mb - b_raw;% Equivalent to E(\hat{\rho} - \rho). Here b_raw is the \rho in alternative-bootstrapped samples
b_adj = b_raw - bias;

% Bias-adjusted adjusted R2
e_adj = py - RHS * b_adj;
R2bsadj = 1 - (sum(e_adj.^2)/sum(dy.^2)) * (T - 1)/(T - size(RHS , 2)); 



Results.b = b_raw; % Unadjusted coefficients
Results.badj = b_adj;% Adjusted coefficients
Results.se = se_raw;   % Unadjusted s.e.

Results.p1alt = p1_alt;   % Bootstrapped p-value under alternative
Results.p1null = p1_null; % Bootstrapped p-value under null

Results.p2alt = p2_alt;   % Bootstrapped p-value under alternative
Results.p2null = p2_null; % Bootstrapped p-value under null

Results.R2adj = R2adj_raw;
Results.R2badj = R2bsadj;
Results.R2nullprc = prctile(R2adj_null , R2prc); % Percentile of adj. R2 under the null

Results.W = XL; %Block length
end